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Abstract 

Random Boolean networks (RBNs) have been a popular model of genetic regulatory net- 
works for more than four decades. However, most RBN studies have been made with random 
topologies, while real regulatory networks have been found to be modular. In this work, we 
extend classical RBNs to define modular RBNs. Statistical experiments and analytical results 
show that modularity has a strong effect on the properties of RBNs. In particular, modu- 
lar RBNs have more attractors and are closer to criticality when chaotic dynamics would be 
expected, compared to classical RBNs. 

1 Introduction 

Random Boolean networks (RBN) have been a popular model of genetic regulatory networks 
(GRNs) [25j [26l US]- Most studies have been made on RBNs with random topologies. Never- 
theless, it has been found that topologies affect considerably the properties of RBNs. For example, 
Aldana studied RBNs with a scale-free topology pQ , discovering important differences with random 

*A preliminary version of this work was presented at the ALife XII conference in Odense, Denmark on August 
20*\ 2010. [36] 
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topologies. In this work, we study the effect of a modular topology in RBNs. We find that mod- 
ularity changes the properties of RBNs. Given the fact that real GRNs are modular [38] EJ 137] 
and most RBN studies have been made over random topologies, it is important to understand the 
differences between random and modular topologies. 

Modularity plays an important role in evolution [42j [14] [50] , since separable functional systems 
are found at all scales of biological systems [48]. Modularity allows for changes to occur within 
modules without propagating to other regions and the combination of modules to explore new 
functions [13] . Thus, the study of modular RBNs is also relevant for understanding the evolution 
of GRNs. 

In the next section, classic RBNs are reviewed, together with their dynamical properties and 
related work. Section [3] presents our model of modular RBNs. Methods and results of statistical 
experiments follow in Section [4] The discussion in Section [5] reflects on the results and provides an 
analytical confirmation. Several future research avenues are mentioned to conclude the paper. 

2 Random Boolean Networks 

Random Boolean Networks (RBNs) [25 , 26, 16 consist of N nodes with a Boolean state, representing 
whether a gene is active ("on" or "one") or inactive ("off" or "zero"). These states are determined 
by the states of K nodes which can be considered as inputs or links towards a node. Because of this, 
RBNs are also known as NK networks or Kauffman models [3 . The states of nodes are decided by 
lookup tables that specify for every 2 K possible combination of input states the future state of the 
node. RBNs are random in the sense that the connectivity (which nodes are inputs of which, see 
Figure [T]) and functionality (lookup tables of each node, see Table [T]) are chosen randomly when a 
network is generated, although these remain fixed as the network is updated each time step. RBNs 
are discrete dynamical networks (DDNs), since they have discrete values, number of states, and 
time [53] . They can also be seen as a generalization of Boolean cellular automata [52] [15] , where 
each node has a different neighborhood and rule. 

RBNs have 2^ possible network states, i.e. all possible combinations of Boolean node states. 
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Figure 1: Topology of an example RBN with N = 5, K = 2. Each node has two inputs that 
determine its state. Since the topology is randomly generated, a node might have several outputs 
or none at all. 



Table 1: Example lookup table to update a node z depending on the state of nodes x and y. Lookup 
tables include all possible combinations of inputs, i.e. 2 K rows. Different nodes will have different 
lookup tables, i.e. Boolean functions. In this example table, the function is the binary XNOR. 
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Transitions between network states determine the state space of the RBN. In classic RBNs, the 
updating is deterministic and synchronous [15]. Since the number of states of the network is finite 
and the dynamics are deterministic, sooner or later a state will be repeated in theory (in practice, 
this can take longer than the age of the universe due to the immense state space). When this 
occurs, the network has reached an attr actor, since the dynamics will remain in that subset of the 
state space. If the attractor consists of only one state, then it is called a point attractor (similar to 
a steady state), whereas an attractor consisting of several states is called a cycle attractor (similar 
to a limit cycle). 

RBNs are dissipative systems, since each state has only one successor, while having the possibility 
of having several predecessor states (many states lead to one state), or no predecessor (a state can 
be reached only from initial conditions, a so called a "Garden of Eden" state). Figure [2] illustrates 
the state transitions of RBNs. 



Figure 2: Example of state transitions. B is a successor state of A and a predecessor of C. States 
can have many predecessors (e.g. B), but only one successor. G is a Garden of Eden state since it 
has no predecessors. The attractor C^D^E^F^C has a period four. [19] 

Note that the topological network (with N nodes, each with a Boolean variable, i.e. one bit) 
is different from the state transition network (with 2^ nodes, each with N bits). One of the 
main avenues in RBN research is involved with studying how the topological network (structure) 
determines the properties of the state transition network (function). 
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2.1 Dynamical Regimes 

RBNs have three dynamical regimes: ordered, chaotic, and critical [53j[16]. Typical dynamics of 
the three regimes can be seen in Figure [3] The ordered regime is characterized by little change, i.e. 
most nodes are static. The chaotic regime is characterized by large changes, i.e. most nodes are 
changing. This implies that RBNs in the ordered regime are robust to perturbations (of states, of 
connectivity, of node functionality). Since most nodes do not change, damage has a low probability 
of spreading through the network. On the contrary, RBNs in the chaotic regime are very fragile: 
since most nodes are changing, damage spreads easily, creating large avalanches that spread through 
the network. The critical regime balances the ordered and chaotic properties: the network is robust 
to damage, but it it not static. This balance has led people to argue that life and computation 
should be within or near the critical regime [30 , 26, 9 , 27 . In the ordered regime, there is robustness, 
but no possibility for dynamics, computation, and exploration of new configurations, i.e. evolution. 
In the chaotic regime, exploration is possible, but the configurations found are fragile, i.e. it is 
difficult to reach persisting patterns (memory). There is recent evidence that real GRNs are in or 
near the critical regime [5]. 

It has been found that the regimes of RBNs depend on several parameters and properties [19] . 
Still, two of the most salient parameters are the connectivity K and the probability p that there is a 
one on the last column of lookup tables. When p = 0.5 there is no probability bias. For p = 0.5, the 
ordered regime is found when K < 2, the chaotic regime when K > 2, and the critical regime when 
K = 2 [12] . The ordered and chaotic regimes are found in distinct phases, while the critical regime 
is found on the phase transition. Derrida and Pomeau found analytically the critical connectivity 

(K c ) = — — ^ , (1) 

V ' 2p(l-p) 1 ' 

This can be explained using the simple method of Luque and Sole [32] : Focussing on a single 
node z, the probability that a damage to it will percolate through the network can be calculated. 



1 This result is for infinite-sized networks. In practice, for finite-sized networks, the precise criticality point may 
be slightly shifted. 
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Figure 3: Dynamics for three RBN with N = 20 and p = 0.5: (A) K = 1 (ordered), (B) K = 2 
(critical), and (C) K = 5 (chaotic). 100 time steps shown from a random initial condition (leftmost 
column, time flows to the right, the state of nodes in time is represented in rows, the RBN state at 
a particular timestep is represented in columns). 
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It can be seen that this probability will increase with as more outputs from a node will increase 
the chances of damage propagation. Focussing on a node j from the outputs of z, then there will be 
a probability p that j = 1. Thus, there will be a probability 1 — p that a damage in i will propagate 
to j. Complement arily, there will be a probability 1 — p that j = 0, with a probability p that a 
damage in i will propagate to j. If there are on average (K) nodes that i can affect, then we can 
expect damage to spread if (K) 2p(l — p) > 1 [32] . which implies chaotic dynamics. This leads to 
the critical connectivity of Derrida and Pomeau [12] , shown in equation [I] 

2.2 Related Work 

We briefly mention particular studies of coupled RBNs, which share similarities with the model of 
modular RBNs presented in the next section. A more detailed comparison can be found elsewhere 

Bastolla and Parisi [6 studied modularity within classical RBNs, i.e. functionally independent 
clusters, but not topological modularity. 

There have been different studies where only two coupled RBNs are considered [22 , 4, 24 . 

There are studies where RBNs are generated in cells of a 2D lattice, similar to a cellular au- 
tomaton, where each RBN is weakly coupled with its von Neumann neighbors [45] [39] [11] . The 
goal is to model intercellular signaling in a tissue. 

Our model is more general than previous models, since there is an arbitrary number of coupled 
networks, and this is not restricted to spatial neighbors. Moreover, it is a natural extension of the 
classic RBN model. 

3 Modular Random Boolean Networks 

We propose a general model of modular random Boolean networks (MRBNs) that extend naturally 
the classic RBN model. A MRBN consists of M modules, each of which is a RBN with N nodes and 
on average (K) (intramodular) inputs per node. Each module has additional (L) (intermodular) 
inputs on average that link random nodes from different modules. These L intermodular connections 
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can be seen as "weak links" [10] between modules. Weak links have been shown to offer stability 
in networks [10]. Figure [4] shows an example MRBN. 




Figure 4: Topology of a example MRBN with iV = 5,M = 4,iT = 2,L=l. Each module is similar 
to the example RBN shown in Figure [I] with one additional input per module (dashed arrows). 
Modules can have several or no outputs. 



Thus, the total number of nodes of a MRBN is given by 



N T = N-M 



(2) 



and the total number of links T is given by 



T = M-((K)-N+(L)) 



(3) 



since each of the M modules has (L) inputs and N nodes with (K) intramodular inputs. 
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The total average number of inputs per node (Kt) is given by 

= £ = & + § (4) 

In the exploration of the space of possible MRBNs, the following measures are useful: 
To study the relationship between number of nodes and modules, the node-to- module ratio /i is 
simply 

To study the relationship between internal (K) and external (L) links, the probability n that a 
link is intramodular is given by 

while the probability A that a link is intermodular is the complement of k\ 

A = 1 - k (7) 

4 Experiments 

The open software laboratory RBNLab [18] was extended to explore the properties of MRBNs. 



RBNLab and its Java source code are available at http://rbn.sourceforge.net 



For all experiments, p — 0.5 and a total number of nodes Nt = 20 was used. Even when this is 
a relatively small size of MRBN, the effects of modularity can be already appreciated. The size of 
networks severely limits the statistical explorations, since each additional node doubles the size of 
the state space S = 2 Nt . 

For each case and each Kt, one thousand networks were generated randomly, exploring one 
thousand randomly chosen initial states for ten thousand steps. This implies at least 10 10 updates 
per MRBN ensemble [28]. 

We performed two sets of experiments: one to explore the statistical properties of different 



9 



families of MRBNs and another to measure the sensitivity to initial conditions of different families 
of MRBNs. In each set of experiments, five cases were studied in two groups. In the first group, 
n (percentage of internal inputs) is explored, while leaving \i (node-to-module ratio) fixed. In the 
second group, \i is varied while (K) = (L). For both sets of experiments, the following five cases 
were considered, varying (Kt) = {1, 2, 3, 4}: 

1. (K) = (L), 

2. (K) = 1, fi -+ 1 

3. (L) = 1, /i 1 

4. TV = 20, M = 1, (L) = , (K) = (K T ) 

5. TV = 1, M = 20, (L) = (if) 

For cases 1, 2 and 3, \i remains fixed (/i = — >• 1) while ft is explored. Case 2 favors intermod- 
ular (external) links, while keeping intramodular links fixed at (K) = 1. Case 3 favors intramodular 
(internal) links, while keeping intermodular links fixed at (L) = 1. Case 1 is intermediate, balancing 
intermodular and intramodular links, restricting (K) = (L). 

Cases 1, 4 and 5 are compared to explore the effect of \i on MRBN properties. Case 4 is equivalent 
to the classical RBN model, since there is only one module (M = 1) and no intermodular links 
((L) = 0). Case 5 is the opposite extreme, where there are N T modules of a single node. In this 
case, intramodular links are self-links, and when (K) > 1 there are in practice fictitious inputs, 
since the dynamic behavior is equivalent to that of (K) = 1. This is because having more than one 
input from the same node is equivalent to having only one input (zero or one), as the same node 
cannot have different states at the same time. Case 5 is different from classical RBNs with only 
self-links since (L) > 0. 

In the following, the variables representing averages such as (K) and (L) will be used without 
the average symbols for simplicity. 
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4.1 Statistical Properties 

For the statistical experiments, the following properties were studied: 

Average Number of Attractors A. This indicates how many distinct sets of states can "at- 
tract" the dynamics of the MRBN. When A > 1 it is considered that the system is multistable 
[44] . There is evidence that in real genetic regulatory networks, attractors correspond to cell 
types [23 , confirming Kauffman's original hypothesis [25] . 

Average Attractor Lengths Le. When Le — 1, there are only point attractors in the network. 
Larger values of Le indicate longer cycle attractors. 

Average Percentage of States in Attractors %SIA. This reflects how much "complexity re- 
duction" is performed by the network, i.e. from all possible 2 Nt states, the percentage of states 
that "capture" the dynamics. Even when complexity reduction is relevant, larger values of 
%SIA indicate a more complex potential functionality of the network, i.e. richer dynamics. 
A large %SIA can be given by large attractor lengths Le and/or a high number of attractors 
A. The more and longer attractors a network has, it can exhibit a richer behavior. 



%S I A =100--^ (8) 

Statistical results often give very high standard deviations a. This is because some networks 
might have a single point attractor, while others might have several cycle attractors. Still, the 
averaged values are informative, showing the effect of different MRBN parameters on the network 
properties and dynamics. Nevertheless, the potential implications of very high standard deviations 
should not be forgotten. For example, in classic RBNs K = N with p = 0.5 implies chaotic 
dynamics. Still, it is possible to construct within this ensemble RBNs with ordered dynamics, e.g. 
having a large number of canalizing functions [2TJ EJ HI] • 
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4.1.1 Number of Attractors 

The A results for cases 1, 2, and 3 are shown as boxplot^Jin Figure |U The primary result is that 
higher values of k yield more attractors. Detailed results (including means and standard deviations) 
and k values can be seen in Tables |AT| |A.2[ and [A3] in Appendix [A] Notice that for K T = 1, case 
2 has the highest k = 1, since the restriction K = 1 implies L = 0. For Kt > 2, case 2 has the 
lowest ft, and also the lowest A. Case 3 has the highest n for Kt > 2, and also the highest A. 




Figure 5: Number of Attractors for different Kt for cases 1, 2 & 3 (K = L, K = 1, and L = 1, 
respectively). Notice logarithmic scale. 



These results suggest that favoring intramodular (if) over intermodular (L) links results in a 
higher number of attractors. This can be explained as follows: if there is a maximum k = 1 L = 0, 

2 A boxplot is a non-parametric representation of a statistical distribution. Each box contains the following 
information: The median (Q2 = xo.50) is represented by the horizontal line inside the box. The lower edge of the 
box represents the lower quartile (Ql = xo.25) and the upper edge represents the upper quart ile (Q3 = xo.75). 
The interquartile range (IQR = xo.75 — a^o.25) is represented by the height of the box. Data which is less than 
Ql — 1.5 • IQR or greater than Q3 + 1.5 • IQR is considered an "outlier" , and is indicated with circles. The "whiskers" 
(horizontal lines connected to the box) show the smallest and largest values that are not outliers. 



12 



i.e. there are no intermodular links. This means that modules are isolated and can be seen as 
independent classic RBNs, with different attractors of different lengths. However, the MRBN will 
consider different combinations of the same modular attractors as different attractors. This can be 
better understood with an example. Let us have a small MRBN M = 2, N = 3, K = 2, L = 0. Since 
modules have no interaction {L = 0), this MRBN is equivalent to two separate classical RBNs. Let 
us assume that the first module has a point attractor: 001 — >• 001 and an attractor of period 2: 
000 — >• 111 — >• 000; and the second module has a point attractor: 000 — >• 000 and an attractor of 
period 3: 100 — >• 010 — >• 001 — » 100. Thus, the combinations of these RBN attractors will yield four 
attractors in the MRBN: 

1. The two point attractors: 001000 -> 001000. 

2. The first point attractor and the period three attractor: 001100 001010 001001 
001100. 

3. The period two attractor and the second point attractor: 000000 — » 111000 — » 000000. 

4. The period two and period three attractors: 000100 -> 111010 -> 000001 111100 
000010 111001 000100. 

Considering that in RBNs A grows algebraically with N [16], MRBNs with several modules M 
and k = 1 will tend to have much more attractors on average than a RBN with the same Nt and 
Kt- This is because the MRBN will have as different attractors all the possible combinations of 
all modular attractors. As intermodular links L are added and n decreases, changes in the states 
of nodes which have L links as outputs might perturb and even destroy attractors. When n is 
minimal, the organization of the MRBN is more similar to a classical RBN, since there are less 
restrictions on where to assign links (in general (for M > 2), there are more possible intermodular 
links (M • (M - 1) • TV 2 ) than possible intramodular links (M • TV 2 )). 

The A results for cases 4, 5, and 1 are compared in Figure [6] It can be seen that /i is also 
relevant for A. More modules (low fi) imply more potential combinations of modular attractors, 
which implies a higher A for MRBNs. Case 4, which is equivalent to a classic RBN has a maximum 
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\i = Nt, i.e. a single module, so there is no possible combination of attractors. Case 4 has the 
lowest A of all five cases. The extreme case 5 has a minimum \i = — — , i.e. N = 1, giving the 
highest A of all five cases. Notice that in case 5 A decreases for Kt > 2. This is because even 
when n is constant in theory, in practice n decreases as Kt increases, given the fact that if K > N 



is equivalent in practice to K = N because of fictitious inputs (TV = 1 in this case, see Table A. 5 
in Appendix [A| . 




Figure 6: Number of Attractors for different Kt for cases 4, 5 & 1 (M = 1, N = 1, and K = L, 
respectively). Notice logarithmic scale. 



4.1.2 Attractor Lengths 

The effect of k on Le seems to be minimal, as the average attractor lengths is very similar for cases 
1, 2, and 3, exponentially increasing with Kt, as it can be seen in Figure [7| 

The effect of \i on Le is seen in Figure [8] For K > 3, classical RBNs (/i = 1) have a higher Le 
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Figure 7: Attractor Lengths for cases 1, 2 & 3 [K = L, K = 1, and L = 1, respectively). Notice 
logarithmic scale. 



than a balanced MRBN (case 1), which have a higher Le than extreme MRBNs with a minimal 
H = —J— (case 5). This can be explained by the fact that in practice attractor lengths grow 
algebraically with N (for deterministic updating schemes, as the one used in this paper) [T7] , For 
the same value of TVt, higher values of \i imply a higher N per module. Having more nodes per 
module allows the possibility of more combinations of states in an attractor, increasing its length. 
It is not so much that a large N favors longer attractors, but a small N restricts their possibility. 

For K < 2, classical RBNs have the lowest Le. Here the modular cases can have longer attractors 
because of the combination of modular attractors is possible with a low L, as explained for A. 
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Figure 8: Attractor Lengths for cases 4, 5 & 1 (M = 1, N = 1, and K = L, respectively). Notice 
logarithmic scale. 



4.1.3 % of States in Attractors 

Since all cases have a constant Nt = 20, the %SIA depends only on A and Le, as shown by 
equation [8] 

For cases where n was varied (1, 2, and 3), it was shown that Le did not vary much depending 
on k. Thus, the results for %SIA shown in Figure [9] are very similar to those of A in Figure [5] 



a higher k, yields a higher %SIA (also seen in the means shown in Tables A.l 1 A.2 1 and |A.3| in 
Appendix [A]) . 



For cases where /i was varied (4, 5, and 1), shown in Figure 10 , the results of %SIA also resemble 
those of A. This is because the differences in the number of attractors (Figure [6| are greater than 
those of attractor lengths (Figure [8]). 
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Figure 9: Percentage of States in Attractors for cases 1, 2 & 3 [K = L, K = 1, and L = 1, 
respectively). Notice logarithmic scale. 

4.2 Sensitivity to Initial Conditions 

One way to characterize the dynamical regime of an ensemble of discrete dynamical systems such 
as MRBNs is by measuring the sensitivity to initial conditions. This is done by measuring how 
small differences in initial states lead to similar or different states: if the system is sensitive to small 
differences, it is considered chaotic. This method is similar to damage spreading [43] or stability 
analysis of dynamical systems [40]. There are also equivalents to Lyapunov exponents in RBNs 
[33] . The rationale is similar in all of these methods: if perturbations do not propagate, then the 
system is in the ordered dynamical phase. If perturbations propagate through the system, then it 
is in the chaotic dynamical phase. The phase transition (critical regime) lies where the size of the 
perturbation remains constant in time. 

To measure statistically the sensitivity to initial conditions of MRBNs, we used the following 
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Figure 10: Percentage of States in Attractors for cases 4, 5 & 1 (M = 1, N = 1, and K = 
respectively). Notice logarithmic scale. 

method [17]: For a randomly generated network, pick a random initial state 5^, and let it run for 
a large number of steps t max (t max = 10000 in our present experiments), to reach a final state Sf. 
Now, apply a random point "mutation" to initial state Si to obtain i.e. do a random bit flip. 
Then, let the network run for t max from S[ to obtain another final state Sf. The difference between 
states can be calculated with the normalized Hamming distance: 

1 n t 

H{A,B) = —Y J \ai-h\ (9) 

i 

If states A and B are equal, then H (A, B) = 0. The maximum H = 1 is given when A is the 
complementary state of 5, i.e. every node with state one in A has a state zero in B and every node 
with state zero in A has a state one in 5, i.e. full anticorrelation. H = 0.5 implies no correlation 
between A and B. The smaller H is, the more similar A and B are. As H increases (up to H = 0.5), 
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it implies that differences between A and B also increase. 

Since there is only one bit difference between Si and S[ and each state has Nt bits: 

Hi = H(Si,Si) = ^- (10) 

Now, to measure the sensitivity to initial conditions, the difference between the final and initial 
Hamming distances AH is used: 



AH = Hf — Hi (11) 

where 

Hf=H(Sf,S' f ) (12) 

A large number of random initial states for a large number of MRBNs are used to calculate an 
average AH for an ensemble. 

If AH < 0, then different initial states converge to the same final state. This is a characteristic 
of the ordered regime, where trajectories in state space tend to converge. If AH > 0, then small 
differences in initial states tend to increase, a characteristic of the chaotic regime, where trajec- 
tories in state space tend to diverge. If AH = 0, them small initial differences are maintained, a 
characteristic of the critical regime, where trajectories in state space neither converge nor diverge 
(in practice, AH w 0). Thus, the average AH can indicate the regime of a MRBN. 

Boxplots of the results for cases 1, 2 and 3 are shown in Figure [TT] Notice that boxplots show 
medians. Means can be compared in Tables pO>l | A. 7\ and |A.8| found in Appendix |A{ 

For Kt = 1, the dynamics are in the ordered regime for all cases, since small differences in initial 
states tend to be reduced, indicated by a neg ative AH = For Kt = 2, Tables [ 



A.6 



A.7 



and 



|A.8| in Appendix [A] show that the average AH is close to zero for all cases, suggesting a critical 
regime. The difference caused by modularity is clearly seen for Kt > 2, i.e. in the chaotic regime. 
The sensitivity to initial conditions is inversely correlated with k. This can be explained as follows: 
the higher the ft, the more "isolated" modules are. Thus, it is more difficult for damage to spread 
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Figure 11: Sensitivity to initial conditions for cases 1, 2 & 3 (K = L, K = 1, and L = 1, 
respectively) . 

between modules, even when average connectivities Kt are high. It can be said that modularity 
in MRBNs brings the dynamics of the chaotic regime closer to criticality. Notice that the phase 
transition (Kt = 2, AH = 0) does not move with k. Rather, the properties of the critical regime 
are expanded into the chaotic regime by higher values of k. 

From Figure [12] it can be seen that \i is also relevant for the sensitivity to initial conditions within 
the chaotic regime (Kt > 2). Case 4 (equivalent to a classical RBN) has the highest sensitivity, 
since damage can equally spread among all nodes in the MRBN. The extreme case 5 has a very 
low sensitivity to initial conditions. On the one hand, it is because of the fictitious links already 
explained. On the other hand, since all nodes have self-connections when K > 1, damage will have 
a lower probability of spreading to other nodes (modules). The intermediate case 1 shows that 
some modularity prevents damage from spreading between modules, bringing the chaotic dynamics 
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closer to the critical regime. 




Figure 12: Sensitivity to initial conditions for cases 4, 5 & 1 (M = 1, N = 1, and K = L, 
respectively) . 



4.2.1 Larger Networks 

To ensure that the results presented above were not an artifact of the small size of the networks 
(Nt = 20) specific experiments were performed to compare cases 3 and 4 with large networks of 
N T 400. For case 3, N = M = 20 and L = 1. For case 4, M = 1, N = N T = 400, if K T , and 
L = 0. For each MRBN family, only one hundred networks were generated and only one hundred 
state pairs were explored for ten thousand steps. These experiments are less statistically significant, 
but they clearly show that the difference in sensitivity to initial conditions is due to modularity 



and not to network size. Results can be observed in Figure [13] It can be seen that the advantage 
of modularity is increased for larger networks. 
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Figure 13: Mean sensitivity to initial conditions for cases 3 & 4 (M = iV, L = 1 and M = 1, L = 0, 
respectively) for large networks [Nt = 400). Error bars indicate standard deviations. 

5 Discussion 

In the previous section, the statistical results showed that MRBNs are more robust than classic 
RBNs for the same Kt values, since modularity reduces the probability of damage to spread. This 
can also be confirmed analytically. 

The method of Luque and Sole [32] can be extended to MRBNs. Instead of focussing on a 
single node with a Boolean state, one can focus on a single module m with internal dynamics. The 
probability that the internal dynamics of m may be perturbed will depend mainly on (L) (as well 
as p). If p = 0.5, then L < 2 implies that on average damage will not propagate between modules, 
even if the internal dynamics of m are chaotic, i.e. K > 2. Still, if K » 2, i.e. m is deep within 
the chaotic phase because of a high number of intramodular links, then the fragility of the network 
will be noticeable in the MRBN (AH » 0) even if damage does not spread outside m. It is clear 
that damage spreading will depend on the value of K as well, e.g. if K is small, damage will have a 
lower probability of propagating within modules, affecting the probability of spread across modules. 

Since damage across the whole MRBN can spread through internal (K) or external (L) links, it 
can be seen that there is the following relationship between critical L and K, extending equation [l] 

{L ^ K ^vh) (13) 
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and for p = 0.5: 



and 



<*W = ^ (16) 



A plot of equation [16] can be seen in Figure [l4j with simulation averages for cases 1, 2, 3, and 
^] Even when the analysis presented above assumes infinitely-sized modules and networks, the 
simulation results with small networks match the theoretical analysis. 



It can be seen that a higher n implies lower L, i.e. lower values in Figure 14 When k = 1, L = 0, 
damage cannot propagate between modules, even for the highest local connectivities (K = N). 
Thus, in principle n can be used to modulate the sensitivity to initial conditions of MRBNs [19] . 

There is a negative correlation between n and Ai7, and a positive correlation between n and A. 
However, the explanations for the effect on sensitivity to initial conditions and on the number of 
attractors seem to be also related to different effects of the topology of MRBNs. Still, it would be 
interesting to study whether it is always the case that RBNs with more attractors on average are 
more robust to damage spread, as it is the case for MRBNs. Alternatively, finding counterexamples 
would be illustrative as well. 



6 Conclusions and Future Work 

We have presented a generalization of random Boolean networks, where modules can be constructed. 
With statistical studies on ensembles of MRBNs and an analytical study, it could be seen that 
different parameters that define MRBNs affect most of the properties of the networks and their 
dynamics. Thus, it can be said that modularity is one way of guiding the self-organization of RBNs 

3 Case 4 is omitted because there is only one module, so the results are not related to how damage propagates 
across modules. Moreover, L = for case 4. 
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0.1 0.2 0.5 1.0 2.0 5.0 10.0 



Figure 14: Theoretical criticality of MRBNs (line) depending on K and L. Above the line there is the 
chaotic phase, while the ordered phase lies below the line. Notice logarithmic scales. Experimental 
averages are shown, except when L = (See Tables A. 6 A.10| ). Different shapes denote different 



cases (see legend), while different colors denote different regimes, obtained from AH (blue for 
ordered, black for critical, red for chaotic; color in online version). Due to finite-size effects (N? = 
20), criticality is considered when —0.0333 < AH < 0.0333. The experimental results that are 
considered as critical can be fitted between the curves L = K _ 2 Q 75 and L = ^q^, shown in gray. 
The deviations from L = are also due to finite sized effects. Results for larger networks are closer 
to L = -| and to AH 0. 

A drawback of the studies of DDN criticality based on sensitivity to initial conditions is that it 
restricts the critical regime to a phase transition. Information-theoretical measures [3TJ |49] might 
offer an alternative to better characterize the critical regime and the effect of modularity and other 
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properties on criticality. 

It was shown here that a high n (high percentage of internal inputs) and a low \i (node-to- 
module ratio) promote criticality in what otherwise would be a chaotic regime. However, is there 
an "optimal" value of n and/or \i for particular systems? It would be also interesting to measure the 
modularity of real GRNs, and measure to what extent the modularity plays a role in their criticality 
[5]. Modularity might help explain why real GRNs tend to have a high average connectivity that 
would set them in the chaotic regime [21], while exhibiting critical behavior [5]. 

The results presented here are encouraging to study modularity at multiple scales, i.e. nested 
modules or hierarchies. RBNs have two scales: nodes and network, with no modularity. MRBNs 
have three scales: nodes, modules, and networks. The MRBN model can be generalized to have 
an arbitrary number of intermediate scales, i.e. modules of modules, with different coupling pref- 
erences. In this way, a multi-scaled modular network could have in principle different dynamical 
regimes at different scales, i.e. damage could propagate at one scale but not at another. A general 
way of defining multiple scales might be with recursive RBNs. 

A scale-free topology has been shown to promote criticality in what otherwise would be ordered 
networks pQ. Scale-free topology and criticality are also present in natural networks [2 [34]. It would 
be interesting to study how modular and scale free topologies could be combined, and whether 
their effects on criticality are cummulative: for abstract RBNs and for real GRNs. For example, 
in our present studies N is constant for all modules, while the module size could have a scale- free 
distribution (few large modules and several small modules). 

The redundancy of nodes [20 and modules [7] has shown to promote robustness in RBNs. 
MRBNs can be general models to study the relationship between modularity, robustness, evolv- 
ability, and criticality. Also, degeneracy can play an important role in robustness [46j |47] and 
evolvability [51], although the role of degeneracy in the criticality of RBNs still remains to be 
explored. 

The potential avenues of research are several. The topics related to modularity and criticality 
are many. We believe that MRBNs can contribute to illuminate these interesting questions. 
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A Tables 



Table A.l: Statistical results for case 1: K = L, \i 1. 
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Table A. 2: Statistical results for case 2: K = 1, \i — » 1. 
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Table A. 3: Statistical results for case 3: L = 1, /i — )> 1. 
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Table A.4: Statistical results for case 4: M = 1, N = 20, L = if = 
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Table A.5: Statistical results for case 5: N = 1, L = K. M = 20. 
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Table A. 6: Sensitivity to initial conditions for case 1: K = L, fi —> 1. 
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Table A. 7: Sensitivity to initial conditions for case 2: K = 1, \i —> 1. 
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Table A. 8: Sensitivity to initial conditions for case 3: L = 1, \i — » 1. 
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Table A. 9: Sensit ivity to initial conditions for case 4: M = 1, N = 20, L = 0, if =if T - 



N 


K 


M 


L 


T 


K T 


AH 






a 


20 


1 


1 





20 


1 


-0.04536 


20 


1 


0.0231 


20 


2 


1 





40 


2 


-0.0111 


20 


1 


0.0930 


20 


3 


1 





60 


3 


0.0964 


20 


1 


0.1806 


20 


4 


1 





80 


4 


0.2381 


20 


1 


0.2076 



Table A. 10: Sensitivity to initial conditions for case 5: N = 1, L — K , M = 20. 
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